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ABSTRACT 

Made-to-measure methods such as the parallel code NMAGIC are powerful tools 
to build galaxy models reproducing observational data. They work by adapting the 
particle weights in an N-body system until the target observables are well matched. 
Here we introduce a moving prior regularization (MPR) method for such particle 
models. It is based on determining from the particles a distribution of priors in phase- 
space, which are updated in parallel with the weight adaptation. This method allows 
one to construct smooth models from noisy data without erasing global phase-space 
gradients. We first apply MPR to a spherical system for which the distribution function 
can in theory be uniquely recovered from idealized data. We show that NMAGIC with 
MPR indeed converges to the true solution with very good accuracy, independent of 
the initial particle model. Compared to the standard weight entropy regularization, 
biases in the anisotropy structure are removed and local fluctuations in the intrinsic 
distribution function are reduced. We then investigate how the uncertainties in the 
inferred dynamical structure increase with less complete and noisier kinematic data, 
and how the dependence on the initial particle model also increases. Finally, we apply 
the MPR technique to the two intermediate-luminosity elliptical galaxies NGC 4697 
and NGC 3379, obtaining smoother dynamical models in luminous and dark matter 
potentials. 
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1 INTRODUCTION 

In galactic dynamics, the modelling of photometric and kine- 
matic observations is of great importance to infer intrin- 
sic properties of galaxies such as their orbital structure, 
total gravitational potential, and phase- space distribution 
functi on (DF). As tau tly summarized bv lSver fc Tremaind 
l|l996l . hereafter different techniques to create made- 

to-measure systems reproducing the observational data have 
been devised, and can be broadly grouped in DF-based, 
moment-based, orbit-based, and particle-based methods. 

DF-based methods fit observations with parametrized 
functions of the integrals of motion or of the action inte- 
grals of orbits. Applications include spherical or integrable 
systems ( e.g. 



[Matthias fc Gerhard|[T999l : lBinnevll2010l V The main advan- 
tage of these methods is that they provide the phase-space 
DF directly, although they generally require assumptions 
on the symmetry of the target galaxy. 

Moment-based methods find solutions of the Jeans 
equations that best reproduce observed quantities 
such as su r face density and veloc i ty dispersion (e.g . 
Young '1980': 'Binney & Mamon' '1982'; 'Binney et al.' 'l990|; 
Magqrrian_&_Binncy ,1994,: Lqkaa 2002: CappcUari 2008 



Williams et"an I2OO9I : Icappellari'et al.l [20091 ). Among the 



1991 



_Gcrhar4 

I995I : iKronawitter et al 



Deionghe|[l98 ^: 'Deion ghe fc de Zeeuwl ll988l: 
Hunter fc d o Zccuw 1 19921 : ICarollo et all 
[2000 ) , axisymmetric models 



(e.g. iHunter fc Qian 



iKuiikenI [T995I : IM, 



agorrian 



1993 



I1995I : 



integrable potentials {e.g. 
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Dehnen fc Gerhard! 1 19941 : 
Merrittlll996l'). and nearly 



Dehnen fc GerhardI 1 199 A 



drawbacks of these methods are the need for assumptions to 
close the system of equations, the lack of any guarantee on 
the positivity of the underlying DF, and the difficulties in 
modelling higher order information such as the line-of-sight 
velocity distribution (LOSVD). 

Orbit-based methods (jSchwarzschildl 1 19791 , Il993l ) 
compute a large library of orbits in a fixed poten- 
tial, and then adjust the weight of each orbit until 
the photometry and kine matics of the target galaxy 
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Ivan den Bosch fc de Zeeuwll2010l ). Schwarzschild modelling 
is widely used, e.g. to infer the masses of black holes at the 
centers of galaxies, but applications are mostly restricted 
to axisymmetric systems. Moreover, the technique requires 
the computation of a large and representative orbit library 
for every new trial potential. 

Particle-based methods for the most part work by 
slowly correcting individual weights of partic les as they are 
evolved in the gravitational potential (jST96l ). until the N- 
body system reproduces the observational data. Kinematic 
and density observables can be used simultaneously in the 
weight correction b y minimizing y'^-deviat ions betwe en dat a 
and particle model (|de Lorenzi et al.ll2007l . hereafter IDLOTI ). 
Among the main strengths of this so-called made-to-measure 
(M2M) technique are its geometric flexibility, the fact that 
the potential can be evolved self-consistently from the par- 
ticles, and that there is no need to store an orbit library. 

The M2M met hod was first applied to the Milky Way's 
bulge and disk in iBissantz et al.l (j2004l ). A version modi- 
fied to model observational data with errors (X^M2M) was 
implemented in the parallel code NMAGIC by IdLOtI . This 
has been used to investigate the dynamics of the outer 
halos of the two intermediat e-luminosity elliptical galax- 
ies NGC 4697 and NGC 3379 (|de Lorenzi et al.ll2008l . BoOgI . 
hereafter IDLO S: DL09)jjind of the massive elliptical galaxy 
NGC 4649 (|Das_et_alJl2QiJ). More re cent imp l ement ations 
of the M2M method can be found in iDehnenI (|2009h . who 
pro posed a different te chnique for the weight adaptation, 
and lLong fc Maol (|2010l ). A related particle method but with 
a different way of adjusting to the observational co nstraints 
is the iterative technique of lRodionov et al] (|2009l '). 

M2M techniques are very promising, but relatively un- 
explored. It is therefore a natural question whether these 
particle methods can actually recover the phase-space DF 
of a target galaxy if the data uniquely specify it. For a given 
set of data, how much does the final particle model depend 
on the initial one? And how is this dependence infiuenced by 
incomplete or noisy data? Furthermore, given that a system 
of A'^ particles is trained to match a much smaller number 
of observational constraints, the problem arises of reducing 
model degeneracies and preventing the method from fitting 
the noise in the data. 

The above issues are related and are connected to 
the concept of regularization. In standard x^M2M prac- 
tice, a weight entropy is used to regularize the particle 
model: through the entropy function all particle weights 
are biased towards a smooth distribution of predefined pri- 
ors, which are specified together with the initial model, 
and thus implicitly contain assumptions about the dy- 
namical structure of the target galaxy. Therefore, unless 
the dynamical structure of the galaxy is known before- 
hand, smoothing with weight entropy makes it difficult 
to construct models with strong phase-space gradients, 
e.g. between near-radial and near-circular orbits. This is 
discussed further in Section (2] A similar effect arises in 
Schwarzschil d models regularized using m aximum-entropy 
constraints (jRichstone fc Tremaind 1988ll, whic h tend to 



isotropize the final DF {e.g. Thomas et al.ll2005l l 



In this paper, we describe a new Moving Prior entropy 
Regularization (MPR) method based on the idea of a distri- 
bution of particle priors, which are computed according to 
phase-space occupation and which evolve together with the 



adaptation of the particle weights. The new method mini- 
mizes the dependence of the solution on the adopted initial 
particle model, and facilitates recovering both a smoother 
and more accurate DF, reducing local fluctuations without 
erasing global phase-space gradients. 

The paper is organized as follows. In Section [2] the ba- 
sics of the x^M2M method are laid out, the main concerns 
related to the traditional regularization are explained, and 
our implementation of MPR is developed. In Section |3] a se- 
ries of spherical target models is constructed for testing the 
M2M method with MPR. Then, in Section |4] and [5] we in- 
vestigate the different roles played by regularization, initial 
particle model, and data quality for recovering the correct 
galaxy model, and we show that the true solution can in- 
deed be recovered from sufficiently good data. Finally, two 
astrophysical applications are presented in Section [S] where 
we reconstruct regularized NMAGIC models for the ellipti- 
cal galaxies NGC 4697 and NGC 3379 in their dark matter 
halos. The paper closes in Section [7] 



2 REGULARIZATION OF PARTICLE MODELS 

In this Section we outline the x^M2M method, and discuss 
some issues related to its standard (weight entropy) regular- 
ization. An alternative method to regularize M2M particle 
models is then presented. For a more detaile d description 
of the M2M t echni que we refer t h e rea der to |ST96| . IDLOTI, 
iDehnenI l|2009l ). and lLong fc Maol (|2010l ). 



2.1 Brief review of x^M2M technique to model 
observational data 

The goal of the x'^M!2M method is to evolve an A''-body 
system of particles orbiting in a potential to make it repro- 
duce the observables of a target galaxy. The potential can 
be either fixed and known a priori, or time- varying and self- 
consistently computed from the particles. 

Each particle is characterized by its phase-space co- 
ordinates Zi — (ri,Vi) and by a weight Wi. The particles 
should be interpreted in a probabilistic sense: they do not 
repre sent single stars but rather p hase-space fiuid elements 
{e.g. iHernauist fc Ostrikeij|l992l ). If M is the total stellar 
mass of the system, then individual particles have masses 

An observable of a target galaxy characterized by a dis- 
tribution function /(z) is defined as 



K,{z)f{z) 



(1) 



where Kj is an appropriate kernel and z — (r, v) are the 
phase-space coordinates. 

Given a set of observables Yj, j = 1, J, including e.g. 
photometry and kinematics, the particle weights Wi of the 
Ai'-body system are evolved until the model observables 



(2) 



agree with the target observables Yj. Here, the kernel in- 
cludes a selection function which ensures that only particles 
with a direct effect on the observable yj contribute to it. 
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Commonly, the model observables are replaced by their 
time-averaged values 

Vj = Vj - ^) dr (3) 

to increase the effective number of particles contributing to 
them, and to reduce temporal fluctuations. 

The task of adapting individual weights of orbiting par- 
ticles until the target and the model observables match is 
accomphshed by solving the set of differential equations re- 
ferred to as "force-of-change" : 

where e is a small positive constant, and the meaning of the 
other variables is clarified below. 

Equation (|4]) maximizes the merit function 

F^-^X^+^^S (5) 
with respect to the particle weights Wi. Here 

is a statistical measurement of the goodness of the fit in 
terms of deviations 

A.W-|(^ (T) 

between target and model observables, taking the error 
(j{Yj) of the target observable into account. 

For the regularization functional, the weight entropy 

JV 

S = - E Uli log(TOi/wi) (8) 

!=1 

is a common choice. S is a measure of the plausibility of the 
model in terms of the smoothness of the weight distribution 
and thus, indirectly, of the resulting DF, and it serves the 
purpose of regularization by pushing the particle weights 
towards some smooth predetermined weights Wi, called pri- 
ors. In typical applications the number of particles is much 
higher than the number of data constraints on the particle 
model; this intrinsic ill-conditioning of the problem trans- 
lates into a large freedom in the weight adaptation, and 
results in models fitting the noise in the data. That is why 
a simple minimization of is not a well-defined procedure 
to determine the model uniquely, and a certain degree of 
regularization is necessary. 

The balance between regularization and fit to observa- 
tional constraints in equation ((5} is controlled by the con- 
stant /i, so that generally models with smaller ^ aim for bet- 
ter fits to the data, but models with larger fj, have smoother 
DFs . In practice, the best choice of ij. is case dependent (se e 
e.g. iGerhard et allll998l : [Thomas et al.|[2005l . lDL08l[DL09l ). 
hinging on the specific properties of the observational data 
to be modelled (error bars, scatter, spatial coverage), the 
phase-space structure of the galaxy, and possibly also the 
adopted initial particle model. 

Finally, note that a likelihood term can be added to 
equation © to account for the constra ints from a sample of 
discrete velocities, as derived in IdlosI . 



2.2 Issues with standard weight entropy 
regularization 

In the framework of the x^M2M method summarized above, 
individual particle weights are slowly adjusted according to 
a compromise between i which pushes them to match the 
target observables, and entropy S, which instead penalises 
against deviations of the weights from the preassigned set 
{wi} of priors; more precisely, from {wi/e} (see equations [5] 
and [8]). 

Even though no rule on the choice of the priors exists, 
they are traditionally set to Wi = wo — 1/N (the "unin- 
formative" or "flat" priors in Bayesian statistics), and the 
same is done for the individual weights of the initial par- 
ticle model. Through the weight adaptation then, the 
standard Global Weight entropy Regularization (hereafter 
GWR) encourages a dynamical structure in the particle 
model which is similar to that of the initial particle sys- 
tem. Of course, this bias is stronger for larger values of /i, 
and wherever the constraining power of the data is smaller, 
e.g. in the outer galactic regions. 

In practice, smoothing the weights globally towards a 
set of preassigned, flat priors through the entropy ([8]) makes 
it difficult to reproduce strong phase-space gradients of the 
target galaxy, e.g. strongly anisotropic velocity distribu- 
tions, unless either the right orbital structure is already in 
place in the initial particle model, i.e. its dynamics is known 
beforehand, or a very small value of ^ is adopted at the 
expense of sm oothne ss of the underlying DF. This w as no- 
ticed both in IdlosI (see their Fig. 10) and IDLOGI , where 
under-smoothed models proved necessary to recover strong 
radial anisotropy in their elliptical galaxy models. 

However, under-smoothed particle models do not rep- 
resent a proper solution. Indeed, sufficient regularization is 
needed not only to prevent the model from fitting the noise 
in the data, but also to oppose fluctuations of the particle 
weights caused by the noise in the data, and to ensure that 
the weight distribution on neighbouring phase-space tori re- 
mains continuous, as intuitively expected for a relaxed stel- 
lar system. 

In what follows we present a new regularization method 
which alleviates the main issues of the standard global 
weight entropy smoothing, and permits smooth M2M par- 
ticle models to be obtained that reproduce the phase-space 
gradients of the target galaxies independently of initial con- 
ditions. 



2.3 Alternative regularization based on moving 
priors 

The logical step forward to ease the issues with the GWR 
is to abandon the idea of constant priors deflned along with 
the initial particle distribution. Instead we will determine 
a smooth distribution of particle priors which follows the 
phase-space structures traced by the weight distribution, as 
the weight distribution evolves to match the observational 
data. Then we will use the weight entropy to bias particle 
weights towards such moving priors. 

This procedure, which we will denote as Moving Prior 
entropy Regularization (MPR) should facilitate a smooth 
DF without erasing larger-scale phase-space gradients. In 
terms of orbits, this means that the new regularization 
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should assign the same prior to particles belonging to the 
same orbit, similar priors to particles on nearby orbits, and 
different priors to particles moving on very different orbits, 
as required in the presence of strong velocity anisotropy. 



2.3.1 Assignment of new local priors 

Therefore, based on Jeans' theorem (e.;;. 
iBinnev fc Tremaind |2008| ). the assignment of new in- 
dividual priors which mirror the underlying evolving DF is 
best based on the integrals of motion, respectively orbits, 
of the particles. This is particularly simple in the spherical 
case, where the integrals of motion are known and can be 
easily found from the particle model. As already pointed 
out bv lDLOgl . the need for regularization is also strongest for 
spherical models, due to their larger number of independent 
orbits with respect to less symmetric systems. The aim of 
the present paper is to show that this method works, and 
how well it works, in the spherical case. A simple axisym- 
metric scheme is shown in Section |B1 and generalizations to 
more complicated geometries are sketched out in Section [T] 
For assigning priors in the spherical case, in practice 
we sort the particles according to their energy E and total 
angular momentum into a rectangular {E, x) grid, where the 
so called circularity integral x = L/Lc is the ratio between 
the actual angular momentum and the angular momentum 
Lc which a circular orbit would have at the given energy E. 
Once particles are binned in a grid of he x rix energy and 
circularity cells, we compute the average weight Wkiik = 
1, . . . , n_B, I = 1, . . . ,nx) contained in each cell, and then we 
assign it as a new prior to all the particles belonging to that 
cell. 

Provided the [E, x) grid correctly resolves the relevant 
phase-space properties of the target, such new priors en- 
sure an orbit-based regularization which acts locally, homog- 
enizing the weights of particles moving on the same and on 
neighbouring orbits, but at the same time tolerates global 
differences among particles on different orbits. 



2.3.2 Smoothing of the grid of particle priors 

We will see that the priors computed in this way can be quite 
noisy. To avoid coarseness in the distribution of priors, and 
so ensure the global smoothness of the underlying model 
DF represented by the {E, x) grid, we implement a two- 
dimensional spline fit of the gridded priors. The technique 
(|Press et al.| [l992') is well kno wn and widely used, also in 
an astrophysical conte xt [e.g. lMerrittlll993l : [Gerhard et al.l 
1 19981 : iDas et al]|2010l ): a thin-plate spline function for the 
new priors on the grid is searched, that minimizes the pe- 
nalized least square function 



prior and its spline value in each {k, l)-ce\\, and 



A^=^e' + A^A(W/), 



(9) 



having defined a function W(E, x) which equals the values 
of the priors on the {E, x) grid. In the equation above, 
measures the deviation between the original value of the 



+ 2 



+ 



dE'^ I \ dEdx j \ dx'^ , ^^^^ 

(10) 

quantifies the complexity of the fitting spline in terms of 
the second derivatives, which are numerically computed via 
finite differences. 

The regularization parameter A determines which of a 
family of splines, ranging from a plane for A — >■ oo to an 
interpolating cubic spline surface for A — )■ 0, is fitted to the 
grid of priors. Obviously, the optimal A is that which resolves 
the relevant structures in the underlying prior distribution, 
but at the same time damps strong and presumably spurious 
variations among nearby priors. 

In principle, A can be calibrated with a sequence of ex- 
periments on the {E, x) grid. However, since particle weights 
evolve in time, and so does the grid of priors, we decided to 
implement t he General Cross Validation technique (GCV, 
IWahba|[l990l ) to determine automatically the optimal value 
of A each time a new grid of priors is computed. GCV is 
based on the principle of sequentially omitting each data 
point, re-fitting the spline, and predicting the value of the 
point from the spline. The technique singles out the optimal 
value of A for this to work best. 



2.3.3 New definition of pseudo-entropy 

The new moving priors substitute the traditional ones in the 
definition Q of the pseudo-entropy, which we slightly mod- 
ify in order to account for the normalization of the weights. 
As already noted, maximizing the standard entropy biases 
the weights towards Wi/e, so that oversmoothing {e.g. for 
high values of /x) causes an undesired global decrease of all 
weights which lea ds to a poor fit of the mass distribution 
(see e.g. Fig. 10 in lDLOSi ). 

In order to avoid such problems, we define a new weight 
entropy 



S - 



log 






(i) 


- 1 









(11) 



for which we can immediately check that (i) maximizing this 
quantity pushes the weights to the actual values of the pri- 
ors, (ii) positive and negative corrections to the weights are 
now a priori equally likely, and (Hi) the whole regularization 
scheme is neutral to mass, so that the only power to alter the 
total mass of the system is left to the data, see Section 13.21 



3 TARGET MODELS AND OBSERVABLES 

In this Section we construct a series of spherical targets to 
be modelled with NMAGIC (Section g] and [5]) in order to 
address two issues, namely (i) testing the ability of the new 
regularization scheme to fit the target data with an intrinsi- 
cally smooth and unbiased particle model, and (ii) exploring 
the extent to which the x'^M2M technique can recover the 
target phase-space structure from a given data set indepen- 
dently of the initial particle model. 

With these aims in mind, we first focus on a prob- 
lem whose solution is theoretically known to be unique. As 
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proved bv lDeionghe fc MerrittJ l|l992l ) in the spherical non- 
rotating case, if the gravitational potential is known and 
complete information on the LOSVD at all radii is available, 
then the underlying DF can be uniquely recovered. There- 
fore, the first target model we design has a known spherical 
potential and is truncated in radius, so that photometric 
and kinematic data can fully constrain it. 

As a second target model, we build an untruncated (in- 
finite) system whose outer regions remain unconstrained by 
data, similar to the case of modelling real galaxies. 

For both target models, we use the 3D luminosity den- 
sity together with the LOSVD along se veral l ong- slits as 
target data in the modelling, similar to IDLOSI and IDL09| . 
For each model, we generate both a set of idealized kine- 
matic data, i.e. a large number of data points with small 
error bars, and a set of more realistic, i.e. sparser and nois- 
ier, data points. We use NMAGIC itself to construct our 
target dynamical equilibrium structures, and to determine 
their observables, as described in more detail in the following 
subsections. 

3.1 Spherical anisotropic Hernquist targets 

Our target models are iHernauistI (|l99d ) spheres with a ra- 
dially anisotr opic orbital st r ucture either of the Osipkov- 
Merritt kind l|Osipkovlll9"79l : lMerrittlll985l . hereafter OM), 
or of the mor e mildly anisotropic, quasi-separable kind 
ijOerhardl Il99lh . Generally speaking, they are isotropic in 
the central regions, and radially anisotropic for radii greater 
than a specified anisotropy radius. 

The potential-density pair for Hernquist models is 



radii r^. The expansion coefficients 



p(r) 



aM 



27rr(r -I- a) 



ip{r) 



GM 
r + a ' 



(12) 



where M is the total mass, G the gravitational constant 
and a the scale length. We set the scale length equal to 1 
kpc, and we use characteristic values of the elliptical galaxy 
NGC 3379 for the total luminosity L = 1.24 x IO^^Lq, the 
stellar mass-to-light ratio T = 5, and the distance D — 9.8 
Mpc. The projected effective radius of our target model is 
i?eff »i 38.3" = 1.82 kpc. 

With respect to the orbital anisotropy, we either fix the 
OM anisotropy radius = 2a, or we use a = 2 and Lq = 
0.3V GMa in the prescription of lGerhardI (|l99lh to generate 
moderately radially anisotropic models (see equations [2.2] 
and [3.14] therein). 

Following the me thod described in 

iDebattista fc SellwoodI (|2000l ). we generate particle 
model realizations of the spherical targets. To construct 
a truncated target, we only retain particles with energies 
lower than iSmax = '/'(''max), with rmax equal to the model 
boundary. 

Finally, we relax the particle models in the fixed Hern- 
quist potential (note that the truncated target is therefore 
not a self-consistent system), and we compute the target 
observables from the final particle model using NMAGIC to 
integrate the particles, as detailed below. 



3.2 Luminosity observables 

We consider as density constraint a spherical harmonics ex- 
pansion of the target luminosity density on a 1-D mesh of 



CLlm,k 



(13) 



are computed from the particle realizations through 
NMAGIC, making use of the clou d-in-cell technique (see e.g. 
IDLOTI . iBinnev fc Tremaind 120081 ') to distribute the weight 
of a particle between nearby grid points. In the definition 
above, L is the total luminosity of the target, Y^^ are the 
standard spherical harmonic functions, and 7^^*^ is the se- 
lection function associated with the cloud-in-cell scheme. 
The radial grid has 60 points, quasi-logarithmically spaced 
between rmin = 0.01" and rmax equal to the model bound- 
ary (for the truncated target) or to 1500" ~ 40-Ro (for the 
infinite target). 

Poissonian error bars, dependent on the number of par- 
ticles in each shell, are assumed for the radial mass, while 
50 Monte-Carlo realizations of the density field of the target 
model allow e rrors t o be assigned to the higher order mass 
moments (see iDLOTh . Because the targets are spherical, all 
model aim,k with Z 7^ 0, m ^ are constrained to be zero 
within these errors, while the aoo.ft are constrained by their 
values for the known target luminosity distribution. 

When comparing the target data with the model observ- 
ables, we compute the latter in the exact same way from the 
particle model. 

By fitting the aim,k coefficients (|13[l . the total luminos- 
ity of the model is adjusted to the target luminosity L. The 
sum of the weights, initially set to X^fli ~ 1> ma^y there- 
fore change if the luminosity of the initial model Linitiai 7^ L. 
In this work, we set Linitiai = L, and we do not adjust the 
mass-to-light ratio, except in Section |6l so that the total 
mass is also constant throughout the evolution. 

3.3 Kinematic observables 

As kinematic target observables, we use the luminosity- 

weighted Gauss-Hermi te mo ments o f the LOSVD 

ijvan der Marel fc Fran:i3 1 19931 : iGerhardI Il993l ) in vari- 
ous slit cells, computed from the particle realizations using 
NMAGIC, through 

^ (14) 



(|DL07l ). Here, Ip is the luminosity in slit cell Cp, dpi se- 
lects only particles belonging to that cell, the dimensionless 
Gauss-Hermite functions are 

u„{u) = (2"+Vn!)"''"i?„(/.)exp(-/.V2) , (15) 
where H„ denote the standard Hermite polynomials, and 
finally 



(«2,i — Vp) /(Jp, 



(16) 



with Vz^i the line-of-sight velocity of particle i, and Vp and (jp 
the best-fitting Gaussian parameters of the target LOSVD 
in the given slit cell. 

To compute Vp, ap, hs, h4 from the particle model, we 
adopt the following procedure. First, we compute the mean 
velocity and rms velocity of particles in each slit cell, and 
use them to estimate the b„^p from the particles through 
equations (|14|l . Next, we use the first order relations 

1 AV ^, 1 Act 
Afti = ; A/12 = F 17 
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-100 100 -100 100 

X (arcsec) x (arcsec) 

Figure 1. Geometry of the two different slit setups. The red 
circle corresponds to iJeff- Left panel: idealized configuration of 
10 slits covering the target; slit cells outside R^f{/2 are square 
(see Section [3.31 1. Right panel: realistic slit configuration, adapted 
from NGC 3379 llDL09h . 



4 CONVERGENCE TO A THEORETICALLY 
UNIQUE SOLUTION 

Here we construct NMAGIC models for the radially 
anisotropic target galaxy model described above. As con- 
straints we use the luminosity density and the idealized 
kinematic data. We determine the amount of regularization, 
investigate whether it is possible to converge to the theoreti- 
cally unique solution, and see how well the target galaxy can 
be reproduced starting from different initial particle models. 

Our NMAGIC models show for this case that if a unique 
inversion of data to recover the underlying target DF exists, 
then it can actually be found via x'^M2M modelling from 
good enough data. The new regularization method proposed 
in this paper significantly improves both the accuracy with 
which the target intrinsic properties are reproduced, and 
the convergence to the right solution independently of initial 
conditions. 



l|van der Marel fc FraM3ll993l : lRix et al.lll997l ') iteratively to 
correct Vp and cTp, until hi and /12 both converge to zero. 
Then, the new target moments bn.p are temporally smoothed 
to reduce fluctuations caused by particle noise, and this leads 
to values of hi and /12 slightly different from zero. Finally, 
the re sulting velocity profile is fitte d by a Gauss-Hermite 
series ()van der Marel fc Franxlll993l ) setting /ii = /12 = 0, 
and the b„^p are recomputed. 

Two different slit configurations are considered, as 
shown in Fig. [T] There, the left panel illustrates a schematic 
view of an idealized slit data setup, which consists of 10 slits 
covering the target and extending as far S-S /"max — 150" ~ 
4JiefE- In order to increase spatial coverage in the outer re- 
gions, and so to decrease the effects of particle noise, slit 
cells outside i?cft/2 are enlarged and made square. More- 
over, Gauss-Hermite coefficients up to ha are considered. To 
assign error bars to the target kinematic data, we compute 
averaged values of the final time-smoothed 6n,p moments in 
the 10 different slits, and then we set the errors equal to y/2 
times the rms deviation of the individual slit cell moments 
from the averag^U- To complete the generation of this slit 
data set, Gaussian random variates with la equal to these 
errors are added to the average moments &„,p. In the follow- 
ing, we refer to this kinematic data set as idealized data. 

For the truncated model, these data are sufficiently 
close to the required "complete" data set that we would 
expect to be able to recover the theoretically unique under- 
lying model with very good accuracy. 

The right panel of Fig. [T] shows instead the 6 slits which 
were used by IDLOGI to model NGC 3379; for this second 
slit configuration, rmax = 100" ~ 3-Rcfi, only t), cr, /is, /i4 
are available, and observational errors for this galaxy are 
adopted. Finally, Gaussian random variates are added to 
the data with la equal to the observational errors. Here- 
after we refer to this latter kinematic data set as realistic 
data. 



^ The factor \/2 in the error bars is necessary because, given this 
generation of kinematic data, the NMAGIC model will have an 
intrinsic particle noise similar to that of the data which it will try 
to match. 



4.1 Modelling procedure and diagnostic quantities 

Starting from an initial particle model, the weights of all 
particles are evolved until the particle system matches the 
target. As initial particle system, we adopt an isotropic 
Hernquist sphere with the same luminosity but scale length 
a — 1.5 kpc. Different velocity distributions are also con- 
sidered, as specified below. During the whole evolution, the 
potential is kept fixed to the target potential. The particles 
are integrated using a leap-frog scheme. 

After a relaxation phase in which the particle system 
is advanced without weight correction, weights are updated 
according to the force-of-change in equation Q, i.e. sub- 
ject to both data constraints and smoothing constraints, for 
~ 10"" correction time steps. We define the model to have 
converged if X^/J averaged over 50 steps is almost constant 
in the last 10* steps, with fiuctuations which are typically of 
order 2%. The particle weights are then constant to a similar 
accuracy with MPR. Finally, the particles are freely evolved 
for another 10* steps without any further weight correction, 
to ensure that the final particle model is well phase-mixed. 
For reference, 10* correction time steps correspond to ~ 42 
circular rotation periods at the target -Rcff • 

When using the new regularization scheme, individual 
priors are not kept constant in time but rather they are 
continuously updated while particle weights are changed to 
match the target observables, as detailed in Section Par- 
ticles are sorted according to their orbital integrals in a grid 
of n_B = 30 and Ux = 10 bins, chosen as a compromise 
between retaining good resolution for the orbit distribution 
and ensuring a sufficient number of particles in all grid cells. 
The average weight contained in each grid cell is computed, 
and then a GCV thin-plate smoothing spline is fitted to the 
distribution of average weights on the grid. The spline value 
in every grid cell is finally assigned as the new prior to all 
particles belonging to the cell. 

A typical outcome of the procedure early in the evolu- 
tion is shown in Fig. [2] where the {E, x) grid of priors is 
plotted before (left) and after (right) smoothing, and simi- 
larly for a horizontal cut (fixed angular momentum) through 
the grid. The cut shows considerable noise before smoothing 
in the central grid cells, where even with a total A'^ ~ 10® 
particles the number of particles per cell for a Hernquist 
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Figure 2. Top: unsmoothed {left) and smoothed {right) grid of 
individual particle priors after ^ 10'^ correction time steps (colour 
bar on the right). Bottom: a cut of the above grid for x = 0.05, 
showing unsmoothed (left) and smoothed {right) particle priors as 
a function of energy. Priors are smoothed among nearby cells with 
the GCV thin-plate smoothing spline described in Section 12.3.21 



model cusp is still small. In the tests presented here, new 
priors are computed frequently in the initial phase and ev- 
ery 10* correction time steps later in the run, which results 
in an efficient regularization at a minimum computational 
cost. When testing the new scheme for very large fi values, 
priors are computed and updated more often. 

The quality of the final particle model of each run is 
assessed through three diagnostic quantities. The first is the 
goodness of the fit to the data in terms of x^/Jj where J 
is the number of data points. Assuming the goodness of fit 
statistics follows a probability distribution function, the 
mean of the distribution is equal to the number of degrees 
of freedom, i.e. the number of constraints (data points plus 
constraints introduced by the merit function) subtracted by 
the number of parameters (model parameters plus fitted 
weights), which are both difficult to quantify. However, if 
the number of degrees of freedom is approximately equal to 
the number of data points, then /J < 1 means that we 
are fitting the data well. 

The second is the level to which the known intrinsic 
kinematics of the target galaxy are recovered by NMAGIC, 
quantified by the rms difference between the intrinsic veloc- 
ity moments of the target galaxy and those of the final parti- 
cle model realization. In the following, the internal kinemat- 
ics (streaming velocity and velocity dispersions) of the parti- 
cle model are computed by binning particles in spherical po- 
lar coordinates, using 21 radial shells quasi-logarithmically 
spaced between rmin = 0.01" and rmax, 12 bins in the az- 
imuthal angle 0, and 21 equally spaced bins in cosO, where 
6 is the polar angle. 

Finally, we determine the degree to which the parti- 
cle model reproduces the known phase-space structure of 
the spherical target. To quantify this we compute the mass- 
weighted relative rms difference between model and target 
weights {wki,m and Wki,t, respectively) in the grid of energy 
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Figure 3. Quality of the final NMAGIC particle models as a func- 
tion of the regularization parameter fi. Top panel: rms deviation 
(%) of the final particle model from the target internal velocity 
moments. Middle panel: rms deviation (%) between the occupa- 
tion of the {E, x) grid of the target and of the final NMAGIC 
model. Bottom panel: goodness of the particle model fit to pho- 
tometric and kinematic data. Crosses refer to models obtained 
with the traditional GWR scheme, dots to models obtained with 
the new MPR. 



and circularity (E, x) used also for the regularization: 

^grid 



E 



■WkLt 



k,l 



(18) 



4.2 Calibrating regularization 

Following the same a pproac h as I Gerhard et al. 
iThomas etal] (120051 ) . IDLOSI . and IDLOQI . we construct 
NMAGIC models for the target galaxy which only differ in 
the adopted regularization scheme and the amount of regu- 
larization, i.e. the value of the parameter jj,. Note that e is 
kept constant between all models. 

The results are summarized in Fig. [3l where the nor- 
malized goodness of fit x^/Jj the mass- weighted rms over 
the {E, x) grid, Ag^id, and finally the rms difference between 
the internal velocity moments of the target and final particle 
model. Akin, are plotted as a function of /i, from unsmoothed 
models (small /i) to oversmoothed models (high fj,). 

We first focus on the NMAGIC particle models obtained 
with the traditional GWR technique (crosses. Fig. [S}. For 
a wide range of values of fi ^ 10*, NMAGIC is able to fit 
the data with x^/J^^- No clear minimum is present in the 
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Figure 4. Intrinsic kinematics of the NMAGIC models obtained 
with the GWR (left panel), and with the new MPR method {right 
panel). From top to bottom: radial velocity dispersion profile, ver- 
tical velocity dispersion profile, and anisotropy parameter. The 
dotted and full black lines show the intrinsic kinematics of the 
initial near-isotropic particle model, and of the truncated target 
galaxy, respectively. Red, green, blue, and light blue lines cor- 
respond to /i = 10"^, 10*, 10''', lO'' adopted in the modelling (see 
Section EjII. 



plotted rms deviations in grid and intrinsic kinematics as a 
function of fi: these quantities stay almost flat for a large 
range of fi, and then rapidly increase for /i>10'*, when the 
increasing amount of smoothing upsets the fit to the data. 
By the time the smoothing becomes effective in damping 
fluctuations in the intrinsic quantities, the bias introduced 
by the global nature of the smoothing has already set in - 
hence no clear optimal value of n is found. For GWR and 
this particular data set, fi = 10* gives a good compromise 
between quality of the data fit and recovery of the target 
properties - but with little smoothing. 

How well the intrinsic kinematics of the target galaxy 
can be recovered is shown in the left panel of Fig. 21 
which compares the known target kinematics with the fi- 
nal NMAGIC models obtained with = 10*, 10^ 10^ 10^ 
As expected, for higher values of fi the internal moments 
remain closer to the initial isotropic moments. 

Fig.[5]shows the level to which the known target DF can 
be recovered by NMAGIC for = 10*, 10®: the distribution 
of total particle weights in the {E, x) grid is plotted for the 
initial particle model, the truncated target, and the models 
obtained with NMAGIC. We denote this by "mass distribu- 
tion function" , or MDF for brevity. Clearly, for fi = 10* the 
main phase-space structures are well recovered, showing that 
NMAGIC is able to fit the data and to approximately find 
the underlying MDF, but the peak on high-E near-radial or- 
bits is underestimated because of the global nature of GWR. 
For the more heavily smoothed case with fi = 10®, this peak 
is completely wiped out. 

We now consider x^M!2M models obtained by fitting 
the same target data with MPR. As can be seen in Fig. [3] 




0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 
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Figure 5. Mass distribution function (MDF) of particle weights 
in the {E,x) grid, for the initial particle model (top left), the 
truncated target galaxy (top right) , the models obtained with the 
standard GWR (bottom) and with the MPR method (middle), for 
different values of /i. The colour scheme reflects the total weight 
contained in each grid cell, where tie = 30, nx = 10. 

(black dots), the new method works very well in reproduc- 
ing the target, and a series of NMAGIC models fitting the 
photometric and kinematics constraints of the galaxy within 
X^/Jjjl can be generated for a wider range of /i values. Of 
these, models obtained for values of fj, <10^ are essentially 
driven by the term alone. However, when regularization 
becomes significant, a minimum is reached both in the rms 
deviation between intrinsic moments of the particle model 
and of the target, and in the rms deviation of their (E, x) 
distributions (top and middle panels of Fig.|31 respectively). 
Remarkably, the minimum is well below that achievable with 
traditional weight entropy smoothing, indicating that the 
phase-space structure and the internal moments of the tar- 
get can be recovered much better using MPR. 

The right panel of Fig. [4] shows how close the in- 
ternal kinematics of the final particle models for n = 
10*, 10*^, 10®, 10'' are to those of the target galaxy. Note that 
the residuals are so small that the trend with seen in Fig. [3] 
cannot be seen; the new scheme allows one to recover the 
target moments almost perfectly. The accuracy with which 
the MDF in {E, x) integral space is reproduced is shown in 
Fig. [S] Visually comparing this plot with the corresponding 
ones for the truncated target and the best model obtained 
using standard GWR, shows that the target is now recovered 
much better. In particular at small energies, i.e. in the outer 
regions, the weight of particles on radial orbits is increased 
while that of particles on circular orbits is decreased more 
effectively with MPR, especially for the preferred fi = 10®. 

We conclude that, for this particular data setup, the 
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Figure 6. Distribution of particle weights for the final opti- 
mally smoothed NMAGIC models obtained with the traditional 
GWR (black histograms, /i = 10*) and with the MPR scheme 
(red histograms, fj, = 10®). Particle weights were initialized to 
wo = 1/750000 ~ 10-6. 



best choice for fj, with MPR is ~ 10 . This value is con- 
siderably larger than the corresponding n of the traditional 
GWR, showing that the new regularization succeeds better 
in reconciling the smoothness of the underlying model with 
orbital anisotropy. 

It is instructive to compare the final distribution of par- 
ticle weights for both regularization schemes. Fig. |6] shows 
that MPR results in a more compact and more structured 
weight distribution, which avoids extended tails of extremely 
increased or decreased weights, while still providing a good 
and less noisy fit the data (see below). A similar comparison 
in the co ntext of Schwarzschild modelling can be found in 
Fig. 17 of lThomas eFaP l|2007h . 

Along with a more compact weight distribution, MPR 
also leads to a smoother particle model. This can be quan- 
tified by computing the rms fluctuations of particle weights 
around the mean value in all the cells of the {E, x) grid for 
the different kinds of regularization, as shown in Fig. [7] 

Not unexpectedly, the fit to the data also looks 
smoother when adopting the new regularization, and the 
larger n value permitted by this scheme opposes an overfil- 
ling of the data points, as can be appreciated for different 
observables in Section |6l 

To summarize, we have tested the x^M2M method with 
a new Moving Prior Regularization scheme for a radially 
anisotropic truncated target model with idealized data, and 
have calibrated the best value of the regularization param- 
eter /i. We have shown that the corresponding NMAGIC 
models match the target data well and recover the MDF for 
this model in its known potential. We have also seen that 
these models are much less sensitive to the value of fi than 
with the traditional weight entropy regularization, which 
can only reproduce the global anisotropy of this model es- 
sentially without smoothing. The new regularization scheme 
allows NMAGIC to recover a particle model that fits the 
data well but is both intrinsically smoother and reproduces 
the properties of the target more accurately. 



Steps 

Figure 7. Rms fluctuations of particle weights around the mean 
value in each {E, x) cell for the optimally smoothed NMAGIC 
models obtained with GWR (black crosses, /i = lO"*) and with 
MPR (red dots, /i = 10^), as a function of the number of 
NMAGIC correction time-steps. The top curve (green trian- 
gles) shows an essentially unsmoothed model (/x = 10^). uiq = 
1/750000 ~ 10~^ is the value of the initial particle weights. The 
grid has ue = 30 times rix = 10 cells, and only those containing 
more than 50 particles are taken into account, to avoid particle 
noise effects in the computation of the residuals. 



4.3 Varying the initial particle model 

The results of the previous subsection already show that 
NMAGIC can recover the orbit distribution of our truncated 
spherical target galaxy from a set of data that specifies it 
essentially uniquely. In this experiment, we used an isotropic 
initial model, so now we investigate the natural question 
whether and how this result is dependent on the choice of 
initial particle model. 

In particular, we consider both the case in which the 
initial particle model has a radially anisotropic OM orbital 
structure (with anisotropy radius ra = 3a, different from 
the target galaxy), and the case in which it is tangentially 
anisotrop i c acco rding to the quasi-separable prescription of 
ICerhardl (|l99ll . with a = 2, Lq = O.sVgM^, c = 0.1). 
For both initial particle models we checked that a minimum 
number of particles on radial and tangential orbits is present 
at each energy. By analogy with the experiment described in 
the previous subsection, the initial particle models are Hern- 
quist spheres with scale-length a = 1.5 kpc, larger than the 
target galaxy (a = 1 kpc). The same setup of the NMAGIC 
run is adopted, together with the optimal /j, values deter- 
mined in Section [4.21 above for the two regularization meth- 
ods. 

With the new MPR scheme, the final NMAGIC models 
obtained for different initial orbital distributions differ re- 
markably little. Table [1] and Figs. 18191 show how well the in- 
trinsic kinematics and phase-space MDF match those of the 
known target galaxy. The MDF of the final particle model 
is very similar to that of the target galaxy, indepentent of 
the choice of initial conditions (Fig. Typical fluctuations 
in the mass-weighted relative rms difference between target 
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Table 1. NMAGIC models for the truncated target galaxy. Differ- 
ent initial particle models are adopted. For the traditional weight 
entropy smoothing pu = 10*; for the new regularization scheme, 
fi = 10^. The goodness of fit I J , Akin and Agi-id are computed 
as described in Section l4.1l In brackets, the same Ajjin and Ag^id 
computed between the initial particle model and the target. 



250 
200 
. 150 
'° 100 
50 
250 
200 
■S.150 
'° 100 
50 
250 
200 
<^ 150 
^ 100 
50 
20 
10 


-10 



0.5 

Cif 

-0.5 



I I I I I ii| 1 — I I I I I ii| 1 — n 



I I 1 1 1 ii| 1 — I I 1 1 1 ii| 1 — h-t 



I I I I I ii| 1 — I I I I I ii| 1 — n 



I I I I I ll| 1 — I I I I I ll| 1 — H 



lJ I i_ 



I I I I lll| 1 — I I I I I ll| 1 — H 




I I I I lll| 1 — I I I I I ll| 1 — h-li 



I I I I lll| 1 — I I I I I ll| 1 — h-li 



0.1 



0.1 



r/R, 



Figure 8. Truncated target with idealized data: recovery of the 
intrinsic kinematics with different initial particle models. From 
top to bottom: radial, azimuthal, and vertical velocity disper- 
sion profiles, mean azimuthal streaming velocity, and anisotropy 
parameter of the NMAGIC models (full lines) for different ini- 
tial conditions (dashed lines). The black line indicates the intrin- 
sic moments of the target galaxy. Blue, red and green colours 
correspond to isotropic, radially anisotropic and tangentially 
anisotropic initial conditions, respectively. GWR was adopted in 
the runs shown in the left panel, while MPR in the runs shown 
in the right panel. 



and model MDF are 12%, while the intrinsic kinematics of 
the target is recovered almost perfectly, as shown in Fig. [S] 
It is instructive to see how a similar result cannot be 
achieved with the traditional weight entropy: an inspection 
of Table [Tl or Fig. [8] and Fig. [9] shows the poorer accuracy 
of the resulting particle models. Especially for the models 
with tangentially anisotropic initial conditions, the smaller 
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Figure 9. Recovery of the MDF of the truncated target (last 
column on the right) for different initial particle models from ide- 
alized kinematic constraints. The first column shows the distribu- 
tion of particle weights in the {E, x) grid for the isotropic, radially 
and tangentially anisotropic initial particle models (from top to 
bottom). The second column corresponds to the final NMAGIC 
particle models obtained with traditional GWR; the third column 
to the models obtained using the new MPR. The colour scheme 
reflects the square root of the total weight contained in each grid 
cell, where we have adopted ue = 30,na; = 10. 



number of particles on radial orbits together with GWR 
makes it more difficult to reproduce the radially anisotropic 
target. 

We conclude that with the new regularization method 
NMAGIC converges to the (theoretically essentially unique) 
solution to a very good level of accuracy, independently of 
the choice of the initial particle model. In this respect, the 
new MPR method is a definite improvement over the tradi- 
tional weight entropy scheme. 



5 EFFECTS OF IMPERFECT DATA 

In astronomical applications, the data constraints are typi- 
cally less stringent than in the idealized case considered so 
far. We therefore now investigate how the results change in 
more realistic circumstances. 

The following tests represent a sequence of problems 
that are increasingly less determined by the data, start- 
ing from the truncated target covered by realistic data 
(Sec. 15. ip . and moving on to infinite stellar systems con- 
strained by data with finite extent. This allows us to isolate 
the different roles played by the quality and completeness 
of data, the initial particle model, and the regularization 
scheme. 

We find that it is still possible to get close to the tar- 
get dynamical structure from different initial particle models 
with the help of the new regularization method, even though 
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the lack and/or poor quality of data introduce degeneracies 
in the models. 



5.1 Truncated target and realistic kinematic 
errors 

First we construct NMAGIC models for a truncated target 
with the realistic kinematic data, with the goal to establish 
how well the target galaxy can then be reproduced from 
different initial particle systems. 

The realistic data have larger error bars and smaller 
data coverage (see Section l3.3p . For these models rmax ~ 
100" is thus smaller than in the previous case. The different 
initial particle models have the same anisotropy structure 
as in Section [4.31 but are adapted to this rmax - hence they 
are more similar to the target. 

We have repeated the analysis described in Section [4.21 
to determine the optimal value of the smoothing parame- 
ter fi when these realistic data are adopted. Results do not 
change much, and suggest that we can keep the values of 
/X = 10* for the GWR and fi = 10"^ for MPR. 

The results of these models are shown in Figs. [TD] (top 
part) and 111! and more quantitatively in Table [T] for both 
regularization methods. The top part of Fig. 1101 shows the 
deviations of the models from the target, comparing the two 
cases in which idealized and realistic constraints are used. 
The three subpanels show the deviations of — cr'^ and — 
a| normalized by the sum of the two velocity dispersions, 
and the deviations of the velocity anisotropy fi, as a function 
of radius. In this figure, the shaded regions correspond to the 
range of deviations obtained from modelling the data with 
the chosen initial particle systems. 

As intuitively might be expected, these deviations in- 
crease (i) moving to larger radii, and (ii) when realistic data 
are considered. The effect of imperfect data on the final mod- 
els is noticeable, in particular closer to the model boundary 
where poorer constraints from slit data worsen the recovery 
of the intrinsic kinematics. 

However, we see that also for realistic data the new 
MPR works well in recovering the internal kinematics of the 
target galaxy independently of the initial particle model, and 
it is superior to the GWR, as deviations are considerably 
smaller. 

The accuracy with which the phase-space mass distri- 
bution function of the target is matched by the NMAGIC 
models is shown in Fig. 1111 where is clear that the traditional 
GWR works less well with these realistic data, especially 
when tangential initial conditions are adopted. 

Because the realistic error bars are larger than those 
used in the experiments in Section |4l the particle weights 
undergo smaller changes until they match target observables 
in a sense [see equations (O and (Q]. For the same reason, 
the final normalized between data and model observables 
turns out to be smaller (see Table [l|. 

To conclude, these experiments show that the new MPR 
method improves both the quality with which the intrinsic 
properties of the target galaxy can be recovered, and the 
independence of the final particle model from the adopted 
initial model. In fact, MPR makes it possible to recover the 
underlying dynamical structure of our radially anisotropic 
target galaxy with good accuracy (A/3 ~ ±0.1 only at the 
outermost point) even when the quality of the data is not 
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Figure 10. Recovery of the intrinsic kinematics for the trun- 
cated and infinite targets (top and bottom figure, respectively) 
with idealized and realistic data (left and right columns, respec- 
tively). The vertical dashed line corresponds to the radial extent 
of the data for the infinite target. The shaded yellow (orange) area 
shows the range of deviations for the MPR, (GWR) method, when 
the specified range of initial conditions is adopted. Full (dashed) 
lines represent the deviations for the final NMAGIC models ob- 
tained with MPR (GWR) starting from different initial particle 
models. Plotted in each panel are, from top to bottom, deviations 
of normalized ct^ — ct^ , cr^ — ct^ , and anisotropy parameter /3 from 
the respective true value of the target; see Sections 15. II and 15. 21 



perfect. As already noticed, MPR allows the use of higher 
H values, thus reducing mass fluctuations and enforcing the 
smoothness of the underlying model without spoiling the fit 
to the data. 



5.2 Finite data for an infinite target 

Real stellar systems are clearly not as sharply truncated in 
radius as the target galaxies studied so far, and their outer 
regions are usually not constrained by the available data. We 
now come to the more realistic case of modelling an infinite 
target galaxy using finite data, to explore the limitations 
that the modelling encounters in this case. 

As target galaxy we consider our usual Hernquist 
sphere with scale length a — 1 kpc, but this time with- 
out truncation. Because of the extreme behaviour of the 
OM radially anisotropic systems at large radii, we choose 
a milder an i sotrop y for our target galaxy, using the models 
of iGerhcirdI (|l99ll ) with specified circularity function [equa- 
tions (2.2) and (3.14) therein], which only depend on a con- 
stant parameter q, set equal to 2. 

Following the same procedure as adopted above, we 
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Figure 11. Recovery of the MDF of the truncated target (last 
column on the right) with realistic data for different initial par- 
ticle models. As Fig. [9] 



model the target starting from different initial particle mod- 
els (isotropic, less radially anisotropic, and more radially 
anisotropic than the target) and using both regularization 
schemes. Idealized data and realistic data are considered in 
turn. 

5.2.1 Infinite target with excellent but radially limited data 

Our modelling of the infinite target with idealized 
but radially limited data confirms previous experiments 
l|Thomas et al.ll2004l . [ra.08) in that the velocity dispersions, 
streaming rotation, and anisotropy parameter of the infinite 
target galaxy can be reproduced reliably only in the regions 
well inside the part of the galax;y covered by the data. 

Quantitative results are reported in Table [51 while the 
bottom left panel of Fig. [TO] shows the range of deviations of 
the final models from the target obtained for different ini- 
tial particle systems. The vertical dashed line corresponds 
to the outermost data point. This panel shows that in the 
inner regions of the galaxy, where the data provide good con- 
straints to the models, the intrinsic properties of the target 
galaxy are well recovered independently of the initial parti- 
cle model, as already found for the truncated target galaxy. 
However, at larger radii, and close to the outermost data 
point, regularization plays a dominant role in the weight 
correction of particles, and in those external regions a bias 
towards the dynamical structure of the initial particle model 
cannot be avoided. 

Nevertheless, our experiments show that the new MPR 
considerably reduces such bias towards the dynamical struc- 
ture of the initial particle model, as can be seen comparing 
the two shaded regions for MPR and GWR. 

If we require |A/3| ^0.1 and compute how far out this 
is achieved for the range of models obtained starting from 
different initial conditions, we find that this radius is 1.4i?cff 
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its fixed potential. Different initial conditions and regularization 
schemes are adopted, as explained in Section 15.21 is the usual 
goodness of fit, normalized by the respective number of observ- 
ables. 

for GWR, while it shifts to 4.3RcS when adopting MPR. 
Considering instead the radius r(A/3 = ±0.2), the standard 
GWR fails at 3.1i?off, while the new method at 8i?eff. 

When, as in this case, a range of dynamical models ob- 
tained from different initial particle models is compatible 
with the data, one could compare and rank models accord- 
ing to the usual goodness-of-fit basis (see e.g. Table [2]), or 
additionally according to a plausibility criterion that, e.g. , 
favours a constant or smooth outer anisotropy profile. 

5.2.2 Infinite target with realistic and finite data 

As a logical final step, we consider the case in which an 
infinite target like the one described above is constrained by 
realistic, rather than idealized, data. 

We model this target starting again from different initial 
particle models, and show the accuracy of the final NMAGIC 
models in the bottom right panel of Fig. IIUI and in Table [21 

The bottom part of Fig. [10] compares the deviations of 
the final models from the target for idealized and realistic 
data. Apparently, the realistic constraints on an infinite tar- 
get galaxy make it really hard for NMAGIC to recover the 
true intrinsic kinematics of the target independently of the 
initial particle model, even though it is still true that the 
new MPR method is superior to the GWR. 

To quantify how well the particle model reproduces the 
intrinsic kinematics of the target galaxy, we can compute 
r(A/? = ±0.1), which is 0.67?efE for the standard GWR, and 
IRcB when adopting the new regularization. Considering in- 
stead the radius r{Al3 = ±0.2), GWR fails at O.QRcH, while 
MPR at 1.4iicft . Here the kinematic data extend to ~ 2J?off. 

Thus, the results previously obtained for the idealized 
data are confirmed: the new regularization scheme provides 
a better reconstruction of the target properties, and is more 
independent on the choice of the initial particle model. How- 
ever, as soon as there is a lack of data to constrain the mod- 
els, regularization becomes the dominant term in the force 
of change acting on particle weights, and the bias towards 
the initial particle model becomes evident. 

The main conclusion from these tests is that the reli- 
ability of our dynamical models is limited to those regions 
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in which good observational data exist, and that the better 
quaUty of the data is reflected in a better recovery of the 
intrinsic properties of the target galaxy. 



6 REGULARIZED PARTICLE MODELS FOR 
NGC 4697 AND NGC 3379 

We now show our new regularization method at work on 
two real galaxies, an d recon stru ct the best-fitting NMAGIC 
models determined in iDLOSl and lDLOgI for the two intermedi- 
ate luminosity elliptical galaxies NGC 4697 and NGC 3379, 

resp ectivel y. 

IDLOSI and IDL09| used NMAGIC to fit spherical and 
axisymmetric models of different inclinations to extensive 
data sets for these galaxies, including photometry, long-slit 
spectroscopy, integral-field data and PNe velocities. Differ- 
ent from the experiments of Sections [4] and (5] particles are 
evolved in a total gravitational potential 

<I>^<P* + <Pd, (19) 

where (j>* is estimated from the A^-particle model for the 
light distribution via a spherical harmonic decomposition 
jSellwoodl [200^ . lDL07l ) assuming a constant mass-to-light 
ratio T, and </)d is a dark matter halo potential with the 
logarithmic parametrization 

2 2 
^o{R,z) = ^lnirf, + R^ + ^). (20) 

Moreover, the mass-to-light ratio is not a fixed parameter, 
but rather it is determined simultaneously with the mod- 
elling of the dynamical structure in the NMAGIC run. 

For both galaxies, the slit data show clear rotation 
along the major axis on the sky. However, the regularization 
scheme as developed in Section 12.3.11 for spherical systems 
discourages any rotation in the particle model, as it biases 
individual weights towards the same prior regardless of the 
sense of rotation of the particles. Thus, in the following tests 
we adopt a modified setup for the grid of priors, binning 
particles according to E and x, and also according to the 
sign of their L2, and assigning individual priors that differ 
between particles with different sense of rotation. For an ax- 
isymmetric potential, this effectively uses the total angular 
momentum L as an approximate third integral, which may 
be expected to be a reasonable fir st approximation unless 
L ~ ~ (|Gerhard fc Sahalll99lf) . 

6.1 The case of NGC 4697 and its halo 

The intermediate luminosity elliptical galaxy NGC 4697 is 
seen almost edge-on. Assuming that the o bserved nuclear 
dust-l ane is settled in the equatorial plane, iDeionghe et al.l 
l|l996i ) derived an inclination i = 78° ±5°, which is consiste nt 
with the bulge-disk decomposition of lScorza et al.l l| 19981 ') if 
the disk component has an intrinsic axis ratio h/R ~ 0.2. 

NGC 4697 has fitted Sersic model index n = 3.5, and 
an effective radius Ras ~ 66" — 3.36 kpc at an assumed 
distance D = 10.5 Mpc. Kinematic data show significant 
major axis rotation reaching ~ 100 km/s at 90". 

NMAGIC axisymmetric particle models assu ming an in- 
clination i = 80° were constructed for NGC 4697 ([dL08) fit- 
ting simultaneously surface brightness photometry, long-slit 



absorption-line kinematics and hundreds of PNe velocities. 
A range of quasi-isothermal halos was found to be consis- 
tent with the observational constraints, and a massive halo 
with circular velocity vo = 250 km /s at 4.3-Rofr, referred to as 
model J in the notation of lDL08l . fits the FN data best. This 
model is characterized by a moderately radially anisotropic 
orbit distribution, with the anisotropy parameter /3 ~ 0.3 at 
the center and increasingly higher in the outer regions. 

These models were constructed using the traditional 
GWR, and the ^ parameter was set to 100 to avoid strong 
biases to the initial conditions. This in turn led to some over- 
fitting of the slit kinematics d ata, es pecially for the higher 
order kinematic moments (see lDLOSi for details). 

Thus we now build a new regularized J model of 
NGC 4697, to see whether a similarly good but smoother 
particle model can be obtained with the help of the new 
MFR. We rerun that exact model with the code NMAGIC 
using the new regularization scheme, and fi = 10^, and us- 
ing constraints from both photometric and kinematic data, 
including PNe. As specified above, we bin particles accord- 
ing to their integrals E,x, and Lz when adopting the new 
regularization method, to allow for the rotation seen in the 
slit data. 

A comparison of the final particle models obtained by 

|DL08| and with this new MPR is shown in Figs. [T^ and [T51 
Fig. [12] shows the projected absorption line kinematics of 
the final par ticle m odels overplotted on the data points. As 
discussed in IDL08| . asymmetries between the left-hand side 
and right-hand side in the profiles do not imply deviations 
from axisymmetry or equilibrium, but rather they are due 
to averages over slightly different slit cells on both sides. As 
expected, these asymmetries decrease when using the new 
MPR, which allows a higher amount of regularizatio n, and 
the model profiles are indeed smoother than for the lDL08l 
model. 

The intrinsic kinematics of the final particle models are 
compared in Fig 1131 the velocity anisotropy increases from 
the center outwards when adopting either GWR or MPR, 
but MPR results in much smoother profiles in the regions 
constrained by data. 

6.2 The case of NGC 3379 and its halo 

In lDLOGl a sequence of spherical and axisymmetric NMAGIC 
models fitting an extensive data set (photometry, long-slit 
and SAURON absorption-line kinematics, FN velocity dis- 
persion data) was constructed to investigate the mass distri- 
bution and orbital structure of the intermediate luminosity 
El galaxy NGC 3379. 

No independent information on the inclination is avail- 
able for this galaxy, and values of i > 40° are consistent 
with the photometry. The effective radius -Refi ~ 47" — 2.23 
kpc at an assumed distance D = 9.8 Mpc. The combined 
kinematic data set shows major axis rotation reaching ~ 50 
km/s at 20", with PNe indicating a further increase to ~ 70 
km /s at 2 20". 

IDLO9I explored a sequence of spherical and axisymmet- 
ric models, together with some triaxial test models. They 
found that their results were insensitive to the adopted ge- 
ometry. Both strongly radially anisotropic models embed- 
ded in massive dark matter halos and nearly isotropic sys- 
tems dominated by the stellar mass are consistent with the 
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Figure 12. Particle model fits to the slit data of NGC 4697 along the major (left) and minor axis (right). The model data points are 
averages over the same slit cells as the target data, and are connected by straight line segments. The black line shows model J of iDLoi. 
while the red line shows the same model obtained using the new MPR. 
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Figure 13. Intern al velocity moments in the equatorial plane for 
model J of iDLoj [black line) and the new regularized J model 
obtained here (red line), for NGC 4697. The vertical dashed lines 
indicate the radial extent of the minor axis slit data, major axis 
slit data, and FN data, from left to right. 



data. I.e., even the extensive data set used in the mod- 
elling was not sufficient to break the mass-anisotropy de- 
generacy (|Binnev fc Mamon"l982'') because of the rapidly 
decreasing velocity dispersion profile for NGC 3379. How- 
ever, an analysis of the quality of the fit and of the likeli- 
hood of the observed PN velocity data for the spherical mod- 
els slightly favoured a range of models centered around the 
radially anisotropic halo C, which was obtained for a quasi- 
isothermal potential in equation (I20|l . with ro/Rctt — 3, 
vq = 130 km/s, and q = 1. 

We now reconstruct that spherical C model with the 
new MPR method. Given the data setup is very similar to 
the one adopted in the tests of Section[5l we set fi — 10®. We 



bin particles according to their integrals E,x, and Lz, be- 
cause of the observed rotation in both the slit and SAURON 
data. For the energy grid we use a linear binning of the func- 
tion exp(i5), which provides a better sampling of the model 
DF in the outer regions for this potential. 

Fig. [14] shows the fit to SAURON data for the final 
NMAGIC models obtained with GWR and MPR. Both par- 
ticle models reproduce the observed rotation with great ac- 
curacy, and the MPR model is clearly smoother than the 
original (symmetrized) data. In particular, notice the ring- 
like structure in the /14 plot. Even though the new model is 
generated using a much higher value of ^, i.e. much stronger 
smoothing, it still does a good job in fitting the observa- 
tional data, with (x^/>.'^)sauron = 0.86 (compared to 0.17 for 
the traditional GWR). 

The intrinsic kinematics of the final NMAGIC models 
(Fig. llSp are similar, but using MPR the kinks in the profiles 
disappear. It can be seen that a strong radial anisotropy is 
required to match the PNe data in this dark matter halo. 



7 DISCUSSION AND CONCLUSIONS 



Building on the work of ST96I. successive investig ations 
(|DL07l : lDL08l:lDL09l : lDehn"enll2009l : lLong fc Madl201oi ) have 
shown the power of the x^M2M modelling technique to learn 
about the dynamics of galaxies. x^M2M methods work by 
adapting the weights of an A^-body particle system until the 
observational data are well matched in a sense, subject to 
additional regularization constraints. These constraints are 
needed to prevent the particle model from acquiring large 
fluctuations because of scatter and noise especially in the 
kinematic data. 

Traditionally, a Global Weight entropy Regularization 
(GWR) is adopted to regularize the underlying particle sys- 
tem. However, through constant fiat priors GWR introduces 
a bias in the particle model which makes it difficult to re- 
produce strong phase-space gradients of the target galaxy, 
e.g. anisotropic velocity distributions, unless its dynamical 
structure is known beforehand. 

In this paper we have described a new Moving Prior 
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Figure 14. Particle model fits to the SAURON integral field kine- 
matic data for NGC 3379. Top row: symmetrized SAURON data. 
Middle row: best-fitting C model (DL09). Bottom row: new reg- 
ularized C model. Mean velocity, velocity dispersion, and higher 
order Gauss-Hermite moments are shown in the panels from left 
to right. 




Figure 15. In ternal velocity moments in the equatorial plane 
for model C of IDLOqI (black) and the new regularized C model 
obtained here (red), for NGC 3379. The vertical dashed line marks 
the last data point. 



Regularization (MPR) method, based on a prior distribu- 
tion for the particles which evolves with the model. Individ- 
ual particle priors are updated along with particle weights 
to keep track of the phase-space structures of the evolving 
weight distribution. The basic idea is to determine the pri- 
ors such that they are similar for particles on neighbouring 
orbits, specified by orbital invariants or integrals of motion 
such as energy and angular momentum in the spherical case. 



The new priors are then used in a weight entropy function 
to ensure a regularization which smoothes locally in phase- 
space, without erasing global phase-space gradients. 

We have then tested this MPR scheme, together with 
the x^M2M modelling technique, using a series of spherical 
target galaxies with both idealized and realistic data. Our 
main conclusions are as follows; 

• For a truncated spherical target galaxy with idealized 
data, for which in theory a unique inversion of the data 
exists, our NMAGIC models with MPR show that the target 
can be recovered accurately, and independent of the initial 
particle model. 

• The new MPR generally improves both the accuracy 
with which the dynamical structure of the target galaxy is 
reproduced, and the convergence to the true solution inde- 
pendent of the initial particle model. Compared to GWR, bi- 
ases in the anisotropy structure are removed, and local fluc- 
tuations in the intrinsic distribution function are reduced. 
Moreover, MPR allows a higher amount of smoothing than 
the weight entropy regularization, while the data are still 
fitted well. 

• Lack or poorer quality of data introduce degeneracies 
in the dynamical modelling results and a dependence on the 
initial particle model, so that the reliability of the models 
is limited to those regions in which good observational data 
exist. Also in this case, the new MPR achieves a better re- 
construction of the target properties and is less dependent 
on the choice of the initial particle model. 

• Using the new MPR, we have reconstructed the best- 
fitting N MAGIC models determined in previous work by 
iDLOsl and lbLOgl for the two elliptical galaxies NGC 4697 and 
NGC 3379 in their dark matter halos. To this goal, we have 
extended the MPR method to the axisymmetric case, using 
the integrals E and Lz and the total angular momentum 
as an approximation to the third integral. The final models 
are intrinsically smoother and provide smoother fits to the 
available data. 

There is clearly room for improving the current version 
of MPR: the method could be generalized to systems of lower 
symmetry using the invariants associated with orbits, e.g. 
the turning points, to assign moving priors in phase-space 
to the particles. Moreover, a cumulative grid-less variant of 
the method could also be implemented. Re-sampling of the 
A'^-body system from ti me to time dur ing and after the ad- 
justment of the weights (lDehnenll2009l ) would enforce equal 
weight for particles orbiting the same torus, but it would 
not take care of smoothing between nearby tori with very 
different weights. 

To conclude, the experiments described in this paper 
show that the moving prior regularization method improves 
the correct and unbiased recovery of the orbit structure of 
the target galaxy from noisy data. A similar regularization 
scheme could also be implemented in Schwarzschild orbit 
superposition models. 
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